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ABSTRACT 

Numerical simulations for the non-linear development of Kelvin-Helmholtz instability 
in two different density layers have been performed with the particle-based method 
(Godunov SPH) developed by Inutsuka (2002). The Godunov SPH can describe the 
Kelvin-Helmholtz instability even with a high density contrast, while the standard 
SPH shows the absence of the instability across a density gradient (Agertz et al. 
2007) . The interaction of a dense blob with a hot ambient medium has been performed 
also. The Godunov SPH describes the formation and evolution of the fingers due 
to the combinations of Ray leigh- Taylor, Richtmyer-Meshkov, and Kelvin-Helmholtz 
instabilities. The blob test result coincides well with the results of the grid-based 
codes. 

An inaccurate handling of a density gradient in the standard SPH has been pointed 
out as the direct reason of the absence of the instabilities. An unphysical force happens 
at the density gradient even in a pressure equilibrium, and repulses particles from the 
initial density discontinuity. Therefore, the initial perturbation damps, and a gap forms 
at the discontinuity. The unphysical force has been studied in terms of the consistency 
of a numerical scheme. Contrary to the standard SPH, the momentum equation of the 
Godunov SPH doesn't use the particle approximation, and has been derived from the 
kernel convolution or a new Lagrangian function. The new Lagrangian function used in 
the Godunov SPH is more analogous to the real Lagrangian function for continuum. 
The momentum equation of the Godunov SPH has much better linear consistency, 
so the unphysical force is greatly reduced compared to the standard SPH in a high 
density contrast. 
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1 INTRODUCTION 



SPH (Smoothed Particle Hydrodynamics, Gingold & Mon- 
aghan 1977 Lucy 19771 is a fully lagrangian and gridless 



method, and has been used widely in the various fields of 



astrophysics (Monaghan 19921, especially, in an irregular 



shaped and/or self-gravitating system. It is because of its 
lagrangian nature and also due to the incorporation of the 



tree-structure (Barnes & Hut 19861. The tree structure is 



very efficient not only in the calculation of the gravity, but 
also in finding neighbours. Therefore, SPH becomes a very 
effective tool in the research of star or galaxy formation. 



However, Agertz et al. (20071 (hereafter A07) showed 



that SPH has a difficulty to describe the Kelvin-Helmholtz 
instability (hereafter KHI) across a density gradient. They 
performed KHI simulations in the two different density lay- 



ers with two standard SPH codes (GADGET2 ( Springel 
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et al. 20021, GASOLINE (Wadsley et al 



grid-based codes (ART 



Kravtsov et al. 



20041 )) and five 
i997f, CHARM 



(Miniati fc Colella 20071, ENZO-PPM (Bryan fc Nor 



man|1997|), ENZO- ZEUS ( |Stone fc Norman|1992[ ), FLASH 
( Fryxell et al.|[2000 |). A complete absence of KHI across a 
density gradient has been observed in the results of the stan- 
dard SPH codes. However, there are nicely rolled vortices in 
the simulations with the grid-based codes even in a high 
density contrast. The standard SPH codes show the vortices 
in the homogeneous density case only. They also performed 
the interaction of a blob and a hot ambient medium with 
a high mach number (the blob test). In the results of the 
grid-based codes, fingers are initiated due to the Rayleigh- 
Taylor and Richtmyer-Meshkov instabilities at the front of 
the compressed blob, and then enhanced by the KHI. Fi- 
nally, the blob is destroyed. However, the standard SPH 
codes show only compression of the blob. They called it "the 
fundamental difference" between the standard SPH and the 
grid-based codes. Their results should be a big problem, be- 
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cause the KHI plays an important role in the various fields 
where SPH has been applied intensively. 

A07 found that there is a strange behaviour of particles 
around the initial contact discontinuity of the two different 
density layers. Particle alignments are observed, and a gap 
forms along the initial contact discontinuity. The gap forma- 
tion or "the peeling" of the particle layers around a density 



gradient has been reported already (Fulk 19941, and the 



modification of the initial particle configuration has been 
suggested as a prescription ( |Monaghan||1987| |Fulk| 11994] ). 
Therefore, A07 performed the same simulation with three 
different initial conditions, such as the lattice, poisson and 
glass. Although the different initial conditions showed dif- 
ferent results to each other, the KHI was still absent in the 
standard SPH simulations. Another possibility as the rea- 
son of the gap formation is the artificial viscosity, because 
the artificial viscosity of the standard SPH has a long his- 



tory of criticism (Watkins et al. 1996 Cha & Whitworth 



2003b I in various aspects. Although the artificial viscosity 



can give a minor change to the results, the KHI is still ab- 
sent irrespective of the artificial viscosity. Furthermore, the 
ENZO-ZEUS code used in A07 employs a von Neumann- 
Richtmyer type artificial viscosity also, but showed the KHI 
in a density gradient. Therefore, the artificial viscosity can- 
not be a reason of the gap formation. Finally, A07 concluded 
that the absence of the KHI and the gap formation are due 
to an inaccurate handling of hydrodynamic force across a 
density gradient. 

The absence of KHI in a density gradient is a serious 
problem to SPH users, so there should be a quick response. 
|Price| ( |2008[ ) suggested an artificial conduction term in the 
energy equation of the standard SPH. Similar to the role of 
the artificial viscosity at the momentum discontinuity, the 
artificial conduction acts on the thermal energy discontinu- 
ity, and changes the pressure profile to a continuous one 
across a density gradient. He showed that the new energy 
equation containing the artificial conduction term can de- 



scribe the KHI across a density gradient. Price (20081 also 
expected the new formulation of Inutsuka ( 2002 hereafter 



102) based on the new Lagrangian function may handle the 
density gradient correctly. 

In this paper, we will revisit the Godunov SPH (here- 
after GSPH) proposed by 102 as a possible solution of the 
inaccurate handling of a density gradient. The same tests in 
A07 have been performed here again with a two-dimensional 
GSPH code. The unphysical force across a density gradient 
is much reduced in the GSPH results, and the KHI and other 
instabilities are observed. Especially, the KHI developing in 
the diagonal direction has been simulated, and a satisfying 
result is obtained. Complicated patterns due to the combi- 
nations of the instabilities develop in the blob test. 

The inaccurate handling of a density gradient in the 
standard SPH has been studied in terms of the consistency of 
a numerical method, and is given in section 2. As a prescrip- 
tion, the momentum equation of GSPH has been revisited 
with the kernel convolution and also the new Lagrangian 
function of a particle system in section 3. The consistency 
of GSPH has been investigated, and a simple test to verify 
the density discontinuity handling has been performed in 
the same section. The KHI simulations in the two different 
density layers and the blob test are in section 4. Finally, the 
summary is given in section 5. 



2 CONSISTENCY OF SPH 

2.1 Stability, consistency and convergence 

Probably, the most important property of a numerical 
scheme is the convergence, because the convergence ad- 
dresses how close a numerical solution is to the actual solu- 
tion. However, it is not easy in general to prove the conver- 
gence of a numerical scheme directly, because the actual so- 
lution is not unveiled in most problems. Therefore, the Lax 



equivalence (or Lax-Richtmyer) theorem (e.g. Gary 1966 
[Ritchmyer fc Morton|[T967l |Despres|[2003l ) is very useful to 
check the convergence of a numerical scheme. According to 
the Lax equivalence theorem, the stability and consistency 
are sufficient conditions of the convergence. 

First of all, the stability can be defined clearly, and has 
been studied intensively in the standard SPH so far JMon- 



aghan|[T989l [Ba lsara 19<^ |Swegle et a"Ll[T995l |Morris| 



199m 



Cha & Whitworth 2003a| at least in the linear regime. The 



standard SPH is conditionally stable, so with the CFL con- 
dition ( Courant fc Friedrichs|1948 1 , the stability of the stan- 
dard SPH is guaranteed 

Secondly, the consistency of a numerical scheme means 
how well the numerical equations of the scheme approximate 



the physical equations (Fulk 1994 LeVeque 20021, and is 



directly related to the analysis of the truncation error. The 
truncation error of a numerical scheme should vanish as the 
time step. At and the grid size. Ax (in SPH, the smoothing 
length, h is comparable to Ax of grid-based codes) approach 
to the infinitesimal value if the scheme has the consistency. 
Therefore, the loss of consistency will lead to low accuracy of 
the numerical scheme. One may concentrate on the consis- 
tency to get the convergence of the standard SPH, because 
the stability is already proved. 

Although the consistency problem of the standard SPH 
is weU known already (e.g. |Fulk|199H |Dilts|1999| ), it wiU be 
reviewed briefiy in the following sections for the convenience 
of readers. Two approximations are needed to get the mo- 
mentum equation of the standard SPH. One is the kernel 
approximation and the other is the particle approximation. 
The consistency of the standard SPH will be examined in 
both of the two approximations. 



2.2 Kernel approximation 

The kernel approximation is given by 



{f){x)= / f{x')Wix-x',h)dx', 



(1) 



where W and {/) are the kernel and kernel-smoothed func- 
tions, respectively. 

In order to check the (order of) consistency of the kernel 
approximation, we will follow the procedure of |Liu fc Liu| 
( 2006) . If a numerical scheme can produce a polynomial of 
up to n*''-order exactly, the numerical scheme is said to 
have the n"'-order consistency. For example, to check the 
0*'*-order consistency of the kernel approximation, put a 
constant function, f{x') = Co into Eq. ([l]), then the kernel- 
smoothed function becomes 



CoW{x — x' , h)dx' = Co- 



(2) 
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Here, the normalisation condition of the kernel function, 

W{x-x',h)dx' = 1 (3) 



is used. The 0*''-order of consistency is easily proved by 
Eq. m. For the f'-order consistency, one may put a linear 
function, Co + C\x' into f{x'), then 

{Co + Cix') Wix - x\ h)dx' = Co + Cix. (4) 



Here, the normalisation condition (Eq. (|3|) and the symme- 
try property of the kernel function, 



/ 



x'W{x — x' , h)dx' — X 



(5) 



are used. 

For a higher order of consistency, we need a higher mo- 
ment of a kernel function. For example, the 2"''-order of 
consistency needs that the 2"^* moment of a kernel function 
should be 



j (x — x'Y W(x — x' ,h)dx' = Q. 



(6) 



However, it is impossible to achieve this with a non-negative 
kernel function. Therefore, the order of consistency of the 
kernel approximation is less than 2 with a non-negative and 
normalised symmetric kernel function. 

For a more complete discussion, we may have to re- 
peat the same analysis above for the first derivative of f{x) 
also, because the hydrodynamic equations contain the first 
derivative of physical quantities. However, it is not necessary 
to check the consistency only, so omitted here. See |Mon-| 
aghan (19821 or Liu et al. (20031 for the further details. 



The consistency depends on the normalisation condi- 
tion of the kernel. Therefore, the kernel approximation loses 
its 0*''-order consistency if the normalisation condition is 
not satisfied, for example, at the edge of a dense cloud in a 
rarefied ambient medium. However, the tests performed in 
A07 and also in Sec. |4] have correct boundary treatments, so 
the incompleteness of a kernel function at the boundary is 
not a critical problem in the tests. 

Although the kernel approximation has the consistency, 
it is not directly used in the equations of SPH. Instead of 
the kernel approximation, the particle approximation is used 
for the derivation of the standard SPH equations, and is 
explained below. 

2.3 Particle approximation 

The particle approximation used in the standard SPH is 
given by 



(7) 



where Wj is W{x — Xj, h). 

We will repeat the same procedure performed in the 
previous section to check the consistency of the particle ap- 
proximation. For the 0"'-order, put a constant Co instead 
of f{xj) of Eq. ([7|, then 



The particle approximation can reproduce the constant 
function when 



E 



1. 



(9) 



Eq. (|9| holds only in an even distribution of particles. There- 
fore, the particle approximation loses its 0*''-order consis- 
tency in an uneven distribution of particles, and eventually 
the standard SPH is unable to converge to the actual solu- 
tion in that situation. 

This problem appears in the equation of motion of the 
standard SPH. With a pressure equilibrium, one of the typ- 
ical motion equation of the standard SPH without the arti- 
ficial viscosity may be written by 



dvi 
~dt 



pt 



+ 



(10) 



where Wij is W{xi — Xj, h), and the physical variables have 
their usual meaning. Although a pressure equilibrium is as- 
sumed, the hydrodynamic acceleration of particle i doesn't 
vanish, so the particle will move. Figure [l] shows this un- 
physical force. The calculated acceleration (red solid line 
with dots) should vanish because the pressure is constant 
across the density discontinuity. However, the acceleration 
shows a repulsion of particles at the discontinuity. This re- 
pulsion damps the initial perturbation, and suppresses the 
KHI. It will make a gap between the two different density 
layers as well. Therefore, one may understand that the oc- 
currence of the unphysical force across a density gradient 
in a pressure equilibrium is due to the loss of the 0*''-order 
consistency. 

The only way to eliminate the unphysical force is to 
make the density term inside the round brackets of Eq. ( 10 1 
an even function. Especially, a uniform density field around 
particle i is the interesting case, and this is why the KHI ap- 
pears in the homogeneous density case (1:1 density contrast 
case) in A07. However, a uniform distribution of particles is 
a special situation, not a general one. 



3 CONSISTENCY OF GODUNOV SPH 
3.1 Kernel convolution 

The inconsistency of the standard SPH is due to Eq. ([9|, 
which appears in the conversion from a continuum (the ker- 
nel approximation) to a particle system (the particle ap- 
proximation). 102 and Dilts (19991 pointed out Eq. ^ as 



a crude assumption, and 102 suggested a density estimation 
at a arbitrary position x, 



With Eq. (Ill, two identities, 

1-E 



p{x[ 



and 



0-E' 



' dx p{x) 



(11) 



(12) 



(13) 



are derived. Instead of Eq. (|9|, Eq. ( [T2| 
cast into the kernel approximation (Eq. 



has been directly 
([T])) to avoid the 
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which is identical to Eq. ([7| if it is evaluated at x = cci. 

In order to check the consistency of GSPH, put a linear 
function, Co + Cix' into f{x') of the right-hand-side of Eq. 
fTil, then 



density 
pressure 

Acc. 
J I L 



-2 3 

u (-x/h) 



Figure 1. The black and blue solid lines are density and pres- 
sure profiles, respectively. The dots on the lines denote the parti- 
cle positions. They implement a density gradient with a pressure 
equilibrium. The particle distance from the centre, x is scaled by 
the smoothing length, h. The red solid line with dots is the ac- 
celeration of particles calculated by the standard SPH. The red 
solid line without dots shows the expected acceleration under a 
pressure equilibrium. A repulsion happens at the density discon- 
tinuity, and will damp the initial perturbation. Therefore, any 
instability across the density gradient may be suppressed in the 
standard SPH. 



inconsistency of the particle approximation. Therefore, the 
final form of the kernel approximation of GSPH becomes 

rCSPH , 



r 



J2 "^J - ^' h)W{x' - Xj,h)dx'. (14) 



p{x' 



The meaning of Eq. ( 14 1 is clear. GSPH considers both host 



and neighbour particles as an extended particle, and uses 
the information of the detailed internal structure of the two 
extended particles. Note that the standard SPH considers 
the host particle as a smoothed one by the contributions of 
its neighbours, bu t th e neighbour particles are still a point. 
jGSPHj-^-j q£ gq_ (2^4 1 reduces to the particle approximated 



function if particle j is considered as a point. UW{x' — Xj, h) 



of Eq. 114 1 is approximated by the delta function, 5(x' — Xj) 
then Eq. ( 14 1 becomes 



E 



(15) 



(Co + Cix') W{x' - x, h)W{x' ~ Xj,h)dx' 



(Co Cix') 
y^^Wix'-x^M 



(16) 



W{x' — X, h)dx' 



= Co + Cix, 



where Eqs. Q, ([5| and (12 1 are used. Eq. (14 1 can repro- 
duce the linear function, so the first order of consistency is 
guaranteed in GSPH. 

Finally, the momentum equation of GSPH is derived 



using Eqs. (Ill - (141, and becomes 
P{x) 



dvi 
~dt 



p2(a;) 



(d,~dj)WiW,dx, 



(17) 



where Wi, di and di are W{x — Xi,h), and respec- 
tively. Contrary to the momentum equation of the standard 



SPH, Eq. (17 1 is expected to converge to the actual solu- 
tion even in a large density gradient. We investigate the be- 



equilibrium to check this. 



haviour of Eq. ( 17 1 at a general density field with a pressure 

{^^-^J)W^W,dx 



dv. 
~di 



= -P 



E' 
E' 



p^{x) 

p{x) p{x) p{x) p{x) 



dx 



-P 
-P 



(km ^ w, dp(x) 

p(x) p^{x) dx 



dx 



(18) 



JL 

dx 



p(x) 



dx 



= 0, 



where Eqs. ( 12 I and ( 13 1 are used. One can see the accelera- 



tion calculated by Eq. (171 vanishes in the pressure equilib- 



rium regardless of the density field, to the degree to which 
the interpolation of the density field used to compute the 
integral in GSPH (see 102) is exact. 



3.2 Perturbation damping test 

In order to observe the particle behaviour in the standard 
SPH and GSPH at a density gradient, a test for the damping 
of a perturbation has been performed. A two-dimensional 
calculation domain, [—L^, Lx] x [—Ly, Ly] has been set. Here, 
Lx = Ly — 7r/2. Particles are located in a lattice by the 
Ax/2 offset initially, and the density contrast is 1 : 2 be- 
tween upper and lower layers. A small displacement of po- 
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Figure 2. The perturbation damping test with the standard 
SPH (right column) and GSPH (right column). Individual snap- 
shot shows the particles position at t = 0.0, 1.3, 2.7 and 4.0tsc 
from the top, respectively. In the standard SPH results, the initial 
perturbation is damped completely in itsc duo to the repulsion 
of particles across the density gradient, although the pressure is 
uniform. In contrast, the perturbation survives in the GSPH re- 
sults. The initial density contrast is 1:2 in the upper and lower 
layers. 



sition, KC^,^!/), 



— — AKsin(Ka;)[exp(— Ky) + exp(Kiy — 2KLy)] 
= — AKCOs(/ta;)[exp(— Kj/) — exp(Ky — 2KLy)], 

~ — An iim{hix)[— exp(Kj/) — exp(— kj/ — 2K,Ly)] 
= Akco8{kx)[— exp(Kj/) + exp(— Ky — 2HiLy)], 



(20) 



is added to the upper and lower layers, respectively. Here, 
X and y are the initial position of the particles (i.e. the lat- 
tice). The amplitude and wavenumber of the position dis- 
placement, A, K are set to O.Olvr and 2, respectively. ^ is 
added to the original position to move the particle to the 
perturbed position. A pressure equilibrium is assumed in 
the whole calculation domain, and the sound speed of the 
upper layer is set to 1. The sound crossing time, tsc of the 
vertical direction in the upper layer becomes 7r/2. Figure [5] 
shows the snapshots of the standard SPH and GSPH results 
at the different times. 

Any movement of the particles is not expected in the 
test, because the pressure is uniform in the whole calculation 



domain. However, the initial contact discontinuity becomes 
flattened in the standard SPH results. The contact discon- 
tinuity of the GSPH results preserves its initial shape well. 
It is clear that the repulsion due to the inconsistency of the 
standard SPH damps the perturbation. The repulsive force 
acts on the normal direction with respect to the density dis- 
continuity, so seems to be a surface tension ( Price|2008 1 . We 
have changed the curvature of the initial perturbation, and 
confirm that the damping depends on the curvature. 



3.3 Lagrangian function 

Another way to derive the equations of SPH is the use of 
a Lagrangian function, (e.g. 'Price & Monaghan 2004 and 



references therein). The Lagrangian function, L of a fluid is 
given by 



/ 



, 1 2 



u I dx. 



(21) 



where u is the specific internal energy. The Lagrangian func- 
tion of the standard SPH is 



(22) 



With this Lagrangian function, the Euler-Lagrange equa- 
tion gives the motion equation of the standard SPH. How- 
ever, the particle approximation is already used in the La- 



grangian function (Eq. (22 1), so the resulting momentum 



equation from the Lagrangian function still has the incon- 
sistency in the uneven particle distribution. 

The relation between the Lagrangian function and the 
exact fluid Lagrangian function is shown in 102. He derived 
the exact Lagrangian function of a particle system, and then 
make an approximated Lagrangian, 



1 .2 

-X, 

2 ' 



u{x)Widx 



(23) 



which has the 2"'^-order accuracy. The new Lagrangian func- 
tion is very similar to that of the standard SPH, but the only 
difference is the specific internal energy term. The specific 
internal energy appears as if smoothed once more than the 
standard SPH, but this form as the second term in the La- 
grangian function is exactly the same as the corresponding 
term in the Lagrangian function for real fiuid (see Eqs. (29) 
and (41) of 102). The momentum equation derived by the 
use of Eq. ( 23 1 is the exactly same as the equation derived 



by the kernel convolution. 

In order to integrate Eq. (171, functional forms of the 
density and pressure are needed. The linear or cubic spline 
interpolation has been used in 102 as the function of the 
density around the particles i and j, but there is a room 
for the further improvement for a more accurate handling 
of the density field. For the determination of the pressure 
and velocity between the particles i and j, a riemann prob- 
lem solver (hereafter RPS) has been used. This is why this 
method is called the "Godunov SPH". As the usual Go- 
dunov grid-based method, any kind of explicit dissipation 
(e.g. artificial viscosity) is not needed by the virtue of the 
RPS. 

Note that the use of an RPS in GSPH has no direct 
relation to either the absence of the KHI or the consistency 



6 Cha et al. 



problem. The unphysical force due to the inconsistent mo- 
mentum equation of the standard SPH has been fixed by the 
new momentum equation of GSPH derived from the kernel 
convolution or the new Lagrangian function. The RPS is 
used for the description of shock waves, because it generates 
a small but sufficient numerical dissipation around shock 
waves. In order to check this point, we have performed the 
KHI simulations with the simplest version of GSPH sug- 



gested by |Cha fc Whitworth (2003aJ . The simplest GSPH 
uses the same momentum equation of the standard SPH, but 
employs an RPS instead of the artificial viscosity. The sim- 
plest GSPH shows also the absence of the KHI in a density 
gradient. 



4 TESTS 

Two kinds of test have been performed. One is the tra- 
ditional KHI simulation in the two layers with a velocity 
shear, and the other is the blob test. All tests have been 
performed with a two-dimensional 2"''-ordei[^ GSPH code 
incorporated with the adiabatic equation of state. The spe- 
cific heat ratio, 7 is set to 5/3 in all simulations. 



Here p„ and pi are the densities of the upper and lower lay- 
ers, respectively, and Vshear is the velocity difference between 
the two layers, tkh is 0.43 in code unit. 

Figure [3] shows the snapshots at different evolution 
times, t = 0.5, 1.0, 1.5 and 2.0tkh- At t = 0.5tkh, the 
initial contact discontinuity is wiggling due to the initial 
perturbation and the velocity shear. There are nicely rolled 
vortices developing around the discontinuity in the later 
snapshots. A distortion of the vortices are observed in the 
snapshot at t — 2.Qtkh, and a mixing layer is expected to 
be formed around the initial contact discontinuity. Finally, 
the mixing layer will stop the KHI. Contrary to the stan- 
dard SPH, GSPH suffers from the unphysical force across 
the density gradient much less than the standard SPH, so 
it can describe the KHI in the different density layers. Note 
that there isn't any kind of additional explicit dissipation, 
such as the artificial viscosity (or artificial conduction) in 
this simulation. 

Figure |4] is the pressure distribution at f = 1.0 and 
2.Qtkh, and shows a less noisy pattern while pressure blips 
are observed across the contact discontinuity in the standard 



SPH result (see figure 6 of |Price| (2008 )). Note that [Price 
(2008| has got a similar pressure map with the artificial 
conduction as well. 



4.1 KHI in the two- layers (p^ : p; = 1 : 2) 

There are two layers with the different density in a pressure 
equilibrium initially. The equilibrium pressure is set to 2.5 
in code unit, and the density ratio between the upper and 
lower layers is set to 1:2. The two layers move to the opposite 
direction to each other with the mach numbers 0.22 and 0.3 
in the upper and lower layers, respectively. The whole cal- 
culation domain is [O, |] x [— |, |] • The size of calculation 
domain is smaller than that of A07 in order to save the cal- 
culation time. The periodic and mirror boundary conditions 
have been implemented in the x and y-directions, respec- 
tively. The total number of particles inside the calculation 
domain is ~ 10^, and the initial configuration of the particle 
distribution is the lattice (A07). 

An initial velocity perturbation in the i/-direction is 
given by 



A„sm(^ — J 



(24) 



where Ao is the amplitude of the perturbation, and set to 
1/40 of the initial velocity shear. Here, A is the wavelength 
of the initial perturbation, and is set to 1/6. Therefore, two 
vortices are expected in the calculation domain. The initial 
perturbation is given only in a thin layer (|j/| < 0.05) around 
the initial contact discontinuity 

With the initial velocity shear and the density contrast, 
the KHI time scale is defined by 



Tkh ~ 



^{pu + pi) 



Vshe 



(25) 



^ The KHI can be triggered with the l*'-order scheme, but 
doesn't develop very well. The details to implement the 2"''-order 
GSPH scheme is very similar to the MUSCL scheme | |van Leer] 
19971, and will be omitted because it is described in 102. 



4.2 KHI in the two-layers [pu : pi = 1 : 10) 

The same KHI simulation presented in the previous section 
has been performed again, but with a different density con- 
trast. The density contrast is much higher than the previous 
simulation, and is set to 1:10. The initial mach numbers are 
set to 0.2 and 0.63 in the upper and lower layers, respec- 
tively. The total number of particles used in this simulation 
is ~ 10^. The initial perturbation is the same as the previous 
simulation. 

Figure [5] shows the results. The two snapshots are at 
t — 1.0 and 1.25tkh, respectively. The earlier stage than 
I.Otkh is very similar to the lower density contrast case 
described in the previous section. However, the vortices are 
not rolled but elongated in the later stage. 

The reason of the elongation is not clear, but we guess 
that it may be due to the poorer resolution of the upper layer 
than the previous simulation ( |Price|[2008| . GSPH (and also 
the standard SPH) is a lagrangian method, so the numerical 
resolution depends on the number density of particles. With 
a similar number of total particles, the higher density con- 
trast between the two layers makes a poorer resolution of 
the lower density layer eventually. Another possible reason 
of the vortex elongation is the initial pressure. We have used 
2.5 as the equilibrium pressure value in this simulation, but 
different choice of the pressure value may change the result. 
However, we'd like to emphasis that the KHI does happen 
in this high density contrast case. 



4.3 KHI in the diagonal direction 

Contrary to grid-based Godunov schemes, in GSPH, all in- 
teractions between the particles i and j reduce to a one- 
dimensional problem on the line joining the two interacting 
particles even in a three-dimensional problem. Therefore, a 
one-dimensional RPS is enough even in a multi-dimensional 
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Figure 3. The KHI simulation in the two different density layers. The initial density contrast between the layers is 1:2, and the initial 
mach numbers of the upper and lower layers are set to 0.22 and 0.3, respectively. The upper layer moves to the right and the lower layer 
moves to the left. The initial contact discontinuity between the two layers begins wiggling due to the initial perturbation, and then the 
nicely rolled vortices develop around the discontinuity. The time of the individual snapshot is normalised by tkh, a^nd shown at the 
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Figure 4. Pressure distribution of the test shown in figure|3|at t = 1.0 and 2.0tx//. These pressure maps are less noisy than the standard 
SPH results (see figure 6 of |Price| | |2008| )). Note that |Pricej |2008| has got a similar result with the artificial conduction term. 



GSPH code. This is an advantage of GSPH than the grid- 
based Godunov schemes, because there is no effective RPS 
in muhi-dimensional situation ( Monaghan|1997 1 . An opera- 
tor sphtting method is essential in the grid-based Godunov 
schemes to describe a multi-dimensional problem with a 
one-dimensional RPS, but any kind of geometrical splitting 
is not needed in GSPH. 

Figure [6] shows the development of the Klfl along the 
diagonal direction. The density contrast is 1:2, and all initial 
conditions are the same as the previous simulation described 
in section |4.1| except the initial particle distribution. The 
initial particle distribution is rotated by 45°. One can see 
the well developed vortices along the diagonal direction in 
the figure. 



4.4 The blob test 

Interactions between dense blobs and strong blast waves are 
an interesting subject in the context of the formation and 



evolution of stars and galaxies (Murray et al. 1993 Klein 



|et al.|p94l [Jones et al.|p96| [Vietri et al.||1997| . If a dense 
blob is exposed to a strong blast wave (e.g. stellar wind or su- 
pernova remnant), the dense blob will be compressed due to 
the blast wave initially, and destroyed finally. The destruc- 
tion of the dense blob is initiated by the Rayleigh-Taylor 
and Richtmyer-Meshkov instabilities ( Inogamov||1999 1, and 
then enhanced by the KHI. However, the instabilities or the 
combinations of instabilities hardly happen in the standard 



SPH due to the unphysical force around the front of the 
compressed blob, so the blob survives for a very long time 
as it is compressed in the hot medium (A07). 

We have performed the blob test with GSPH. The cal- 
culation domain of the blob test is [—2, 30] x [—6, 6] in code 
unit. A dense blob is at the origin initially, and is surrounded 
by the hot ambient medium moving in the x-direction. The 
radius of the blob is 1, and the density ratio between the 
ambient medium and the blob, x is set to 10. The initial 
mach number of the ambient medium is 5. The numbers of 
particles to implement the blob and the ambient medium 
are 7688 and 93139, respectively. The initial configuration 
of the particle distribution is the glass (A07). The sound 
speed and the density of the ambient medium are set to 1. 



With this initial condition, the cloud crushing time (Klein 
et al. 1 1994 1, Tcc is determined by 



(26) 



where rh, and Va are the radius of the blob and the veloc- 



ity of the ambient medium, respectively. Jones et al. ( 1996 1 



defined the "bullet crushing time", but the only difference 
between Tcc and the bullet crushing time is a numerical fac- 
tor (= 2), so we have used Tcc as the time unit in the blob 
test. Finally, the KHI time scale (A07) of the blob test is 
defined by 



TKHI,blob ~ 1.6 X 2Tcc- 



(27) 
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Figure 5. The KHI simulation with a much higher density contrast. The lower layer is 10 times denser than the upper layer. The vortices 
develop around the discontinuity, and then elongate in the long— term evolution. The elongation may be due to the low resolution of the 
upper layer. 



Tec of this blob test is 0.63. 

Figure [7] is the result of the blob test. The early evolu- 
tion stage of the interaction is compression. The front of the 
blob is compressed due to the ram pressure of the ambient 
medium. There is an evaporation behind the blob. A bow 
shock forms around the compressed blob, and then three 
fingers develop due to the Rayleigh-Taylor and Richtmyer- 
Meshkov instabilities. The fingers are enhanced by the KHI, 
so the mushroom pattern develops at the head of the fingers 
(e.g. Yabe et al.|1991 ). The result of the blob test performed 
by GSPH is similar to the results of the grid-based code (e.g. 
Klein et al.||1994[ ). 



5 SUMMARY 

The standard SPH does not accurately describe pressure 
gradient in the location with a large density gradient, so it 
shows the absence of the KHI in that situation. This is due to 
the inaccurate force calculation across the density gradient. 
There is an unphysical force across the density gradient, and 
this unphysical force pushes the particles away from the the 
initial discontinuity to make a gap and to damp the initial 
perturbation. Therefore, the development of any instability 
is suppressed at the density gradient. 

The inaccurate force calculation is due to the incon- 
sistency of the standard SPH. The particle approximation 



used in the derivation of the motion equation of the standard 
SPH loses the 0*''-order consistency if particles are unevenly 
distributed. One may use the Lagrangian function for the 
derivation of the momentum equation, but the Lagrangian 
function of the standard SPH uses the particle approxima- 
tion already, so the resulting momentum equation shows still 
the unphysical force across a density gradient. 

In order to solve the consistency problem of the stan- 
dard SPH, we have revisited the new formulation of 102, 
called GSPH. With the kernel convolution, new momentum 
equation is derived. We have proved that the momentum 
equation of GSPH has linear consistency up to the accuracy 
with which the kernel convolution integral can be calculated, 
leading to a much reduced unphysical force across a density 
gradient in the pressure equilibrium. The same momentum 
equation can be derived using the new Lagrangian function 
(102). It is very similar to the Lagrangian function of the 
standard SPH, but is more accurate to the real fiuid La- 
grangian function. 

We have explained the geometrical meaning of the mo- 
tion equation of GSPH. It considers the host and neighbour 
particles as an extended body, and uses the detailed informa- 
tion of the extended bodies. In the standard SPH, the host 
particle is smoothed by the contributions of neighbours, but 
the neighbours are considered as a point. 

Two kinds of test have been performed to show the 
performance of GSPH. One is the traditional KHI test in 
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Figure 6. The KHI developing in the diagonal direction. This is essentially the same test presented in figure [3] but the initial particle 
distribution is rotated by 45°. Contrary to grid— based Godunov methods, GSPH can describe a multi— dimensional problem with a 
one-dimensional riemann problem solver, so any kind of operator splitting is not needed. 



the two layers, and the other is the blob test. In the two 
layer test, GSPH showed the development of the KHI even 
with the very high density contrast. The KHI developing 
along the diagonal direction has been performed also, and a 
satisfying result has been obtained. In the blob test, GSPH 
can describe the formation and evolution of the fingers due 
to the instabilities and the combinations of instabilities in 
front of the compressed blob. The blob test result of GSPH 
coincides with the results of the grid-based codes. 

In the standard SPH, not only the momentum equation, 
but also the energy equation is inconsistent in the uneven 
particle distribution. We are investigating the influence of 
the inconsistency on the energy equation, and it is left for a 
following work. 
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Figure 7. Result of the blob test. Each snapshot is at t = 1, 3, 4, 5, 6 and Ttcc from the top-left, and shows the different stage of 
the blob evolution inside a hot moving ambient medium. We have taken the square root of the column density to generate the surface 
density color map. The initial stage is the compression of the blob due to the ram pressure of the ambient medium. One can see the 
bow shock formed in front of the compressed blob. The ambient medium enters into the bow shock is decelerated, and then eventually 
the Rayleigh-Taylor and Richtmyer— Meshkov instabilities are triggered. Three fingers begin to develop in front of the compressed blob 
in the third snapshot, and then the fingers are enhanced by the combinations of the instabilities. A mushroom pattern develops at the 
head of the fingers. This result coincides with the results of the grid— based codes well. Note that this is a two dimensional simulation, 
while the blob test of A07 has been performed in three dimensions. 
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